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The method for analytic evaluation of four-particle integrals, proposed by Fromm and Hill, is 
generalized to include complex exponential parameters. An original procedure of numerical branch 
tracking for multiple valued functions is developed. It allows high precision variational solution of 
the Coulomb four-body problem in a basis of exponential-trigonometric functions of interparticle 
separations. Numerical results demonstrate high efficiency and versatility of the new method. 
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I. INTRODUCTION 

The problem of four particles with the Coulomb in- 
teraction plays an important role in atomic and nuclear 
physics. It forms a link between the three-body prob- 
lem, which can be solved numerically with very high 
precision, and many-body problems, solutions of which 
are very approximate. Thus, profound studies of various 
four-particle systems can provide valuable insights into 
physics of systems with greater numbers of particles. 

In addition to the methodological interest, the four- 
body problem has unquestionable practical significance. 
Positronium beams are extensively used in positronium- 
atom scattering experiments, but the positronium 
molecule, e+e^e+e^, has not been observed experimen- 
tally yet. All existing knowledge of its properties is based 
on numerical studies [Q. Molecules and ions including 
/i-meson have attracted much attention traditionally in 
connection with the problem of the muon catalyzed fu- 
sion. Calculations suggest that muonic molecules like 
p^fi^p^fi^ have higher nuclear reaction rates than the 
corresponding three-particle ions. These examples show 
that high-precision numerical solution of the four-body 
problem is essential for proper understanding of various 
physical phenomena. 

The majority of four-particle systems are nonadiabatic, 
and cannot be treated within the adiabatic approxima- 
tion. The only practical way to calculate their energy and 
properties is to use the variational approach, taking into 
account the correlated motion of all the particles. Ba- 
sis functions of the Gaussian type, depending on six in- 
terparticle separations and several nonlinear parameters, 
have been extensively used for such calculations ||. |j . 
An important advantage of the Gaussian functions is that 
all integrals can be easily evaluated. The nonlinear pa- 
rameters are optimized stochastically ||^ : at each step of 
basis expansion, many functions with randomly gener- 
ated parameters are examined, and the function, giving 
the largest decrease in energy, is added to the basis. 



However, unlike real wavefunctions, the Gaussian func- 
tions do not decay exponentially, and do not satisfy the 
cusp condition. From this point of view, they are rather 
unphysical. As a result, convergence of the variational 
procedure is very slow, and many hundreds of basis func- 
tions must be used. A recent calculation of the positro- 
nium molecule ||] involved 1600 Gaussian functions. It 
was suggested that further expansion of the basis was not 
practical because of increasing computation time and low 
probability of finding good parameters. Thus, more effi- 
cient basis functions are clearly required. 

A method for analytic evaluation of four-particle inte- 
grals, proposed by Fromm and Hill opened up possi- 
bilities of variational calculation of four-particle systems 
in a basis of exponential functions of interparticle sepa- 
rations. This method reduces computation of integrals, 
needed to determine matrix elements of a four-particle 
Hamiltonian, to evaluation of the dilogarithmic function 
Q of various arguments. Application of this method, 
however, is a very difficult problem. Because the diloga- 
rithm is a multiple valued function, the entire algorithm 
cannot be used without an effective procedure of branch 
and singularity tracking. 

This problem was initially solved by the authors for 
the case of real exponential parameters. The first calcu- 
lations of the positronium molecule |^ , and several mesic 
molecules Q in the exponential basis, depending on all 
six interparticle separations, have demonstrated high ef- 
ficiency and great potential of this method. To the best 
of our knowledge, nobody else has done this yet 

Because one exponential function is as effective as eight 
Gaussians, a size of the basis can be reduced significantly. 
However, an amount of time, needed to compute one ma- 
trix element, is much larger than for the Gaussian basis. 
Thus, optimization of nonlinear parameters is the main 
difhculty. Deterministic optimization (gradient descent) 
gives excellent results for a relatively small number of ex- 
ponential basis functions. Stochastic optimization (trial 
and error), used to expand the basis further, is ineffi- 
cient due to a dramatic increase in computation time. 
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This fact suggests that a possible alternative to an enor- 
mously large Gaussian basis is a relatively short basis of 
the most efficient and versatile functions with carefully 
optimized parameters. 

A natural generalization of the exponential basis is 
the exponential-trigonometric basis, obtained by replac- 
ing real exponential parameters with complex ones 
The exponential-trigonometric functions have been suc- 
cessfully employed in variational calculations of three- 
particle adiabatic systems jl^. They are much more ef- 
ficient, than the ordinary exponentials for two reasons. 
First, they contain twice as many nonlinear parameters, 
thus allowing better approximation of the wavefunction. 
Second, they exhibit nonmonotonic dependence on inter- 
particle separations, being able to imitate sharp peaks 
in wavefunctions of adiabatic systems. The computation 
time increases only insignificantly in comparison with the 
case of real exponential parameters. 

In order to use the exponential-trigonometric basis in 
the four-body problem, one has to evaluate the four- 
particle integrals with complex parameters. The prob- 
lem of branch tracking in a general complex case is 
formidable. Every branch change for every multiple val- 
ued function has to be taken into account if correct val- 
ues of the integrals are to be obtained. An original 
(and, inevitably, very nontrivial) procedure of numeri- 
cal branch tracking has been developed by the authors. 
The first variational calculations of four-particle systems 
in the exponential-trigonometric basis proved extremely 
promising [ pl| . They showed that one exponential- 
trigonometric function can replace seven exponential 
functions in calculation of e^e~e~^e~ , and several dozen 
exponentials in studies of adiabatic systems . There- 
fore, it presents a real alternative to both the exponential 
and Gaussian basis functions. 

Even though results of the calculations involving 
the exponential and exponential-trigonometric functions 
have been published ^, ^, ^ , details of the new method 
have not been reported yet. The purpose of the present 
paper is to fill this gap. We present a description of our 
algorithm that will enable a reader to implement it as a 
computer program. 

The paper is organized as follows. Section II. A dis- 
cusses what integrals are needed to compute matrix ele- 
ments of a four-particle Hamiltonian, and how the num- 
ber of them can be reduced. In Sec. II. B, principles of 
the original method by Fromni and Hill are outlined. 
Sec. II. C provides information about multiple valued 
functions used in the analysis. In Sec. II. D, a simplified 
procedure of branch tracking in the case of real parame- 
ters is described. Sec. II. E gives a detailed exposition of 
the method of branch tracking in the most general case, 
when all the parameters are complex. In Sec. II. F, a prac- 
tical implementation of the branch tracking algorithm is 
described. The last Section presents our conclusions. 



II. DESCRIPTION OF THE METHOD 

A. Matrix elements of four-particle Hamiltonian 

Let us consider a Hamiltonian of a four-particle system 
with the Coulomb interactions: 



H 
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3 = 1 



j<k 



(1) 



Here nij and q^, j = 1...4, are masses and charges, and 
rjk — \rj — Tfcl are interparticle separations. Our purpose 
is to evaluate matrix elements of H with exponential ba- 
sis functions 

4 4 

*6 = cxp(- ^ bjkVjk), $c = cxp(- ^ Cjkrjk)- (2) 

'j<k j<k 

These functions depend on complex parameters {bjk} 
and {cjk}- In what follows, the notation {xjk} will al- 
ways refer to six quantities, Xjk, with j,k — 1...4 and 
j < k, i.e. a;i2, Xia, a;i4, a;23, a;24, 2^34 7 assuming that 

Xjk — -^kj . 

In order to compute matrix elements of the operator 
of kinetic energy in Eq. (1), one has to evaluate the fol- 
lowing quantities: 



($h I cos 9 



jki 



y->2 y-»2 



kl 



2rikr 



jk"jl 



l*c), (3) 



where j ^ k, k ^ I, j I. The integrands in the last 
formula display linear and even quadratic dependence on 
certain interparticle separations. Therefore, in order to 
obtain matrix elements of the Hamiltonian, Eq. (1), one 
has to calculate a total of 43 integrals: one overlap in- 
tegral, six integrals of the Coulomb interactions, and 36 
integrals, given by Eq. (3). 

It turns out, however, that it is possible to avoid com- 
putation of the integrals in Eq. (3). It has been shown by 
one of the authors that the matrix elements of the above 
Hamiltonian can be expressed in terms of the overlap in- 
tegral and six Coulomb integrals only . Thus, one can 
write: 



($ J I =H,-H2-H. 



(4) 



The individual terms in Eq. (4) are given by the following 
expressions |l2[: 



j<k 



{m.j H- mfc ) 
2mjnik 



-ajk + QjQk 



($J — I ; (5) 
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J=l k<l 



{ajkSjk + ajiSji — ajnSjn) 
2mj ajkaji 



djkdji . (7) 



In these formulas, $a is a new function with parame- 
ters {fljfc}, defined as a^k = {bjk + Cjfe)/2: 



$a = exp(- E o.jkrjk) 
j<k 



(8) 



The parameters {djfe} are defined as djk ~ (cjk — bjk)/2, 
and the quantities {sjk} are given by: 

Sjk = ($a I — I $a) - ajfc($a | $a) • (9) 

The additional index ri in Eq. (7) is fixed by a condition 
n ^ j, k, I. 

Therefore, only seven integrals - the overlap integral 
and six Coulomb integrals, calculated with the function 
$a, - are needed to determine the matrix elements of the 
Hamiltonian, Eq. (1). The above formulas are valid for 
both real and complex parameters. They are indispens- 
able for any application of this method. 



B. Evaluation of four- particle integrals 

In this paper, we generalize the method of analytic 
evaluation of four-particle integrals, proposed by Fromm 
and Hill Q, to include complex exponential parameters. 
First, we would like to recall basic ideas of this method. 
The following family of integrals is considered: 

•^({"jfe}){ajfc}) = \(S\"^7k ^)exp(-E"jfe^jfe) ■ 



j<k 



j<k 



(10) 

Here, {ajk} denotes a set of six exponential parameters, 
"12, ai3, "14, ce23, Q;24, ^34, and {ujk} is the corre- 
sponding set of non-negative integers. The integrand de- 
pends on six interparticle separations, {rjk}- The inte- 
gration is performed over 9-dimensional space of relative 
coordinates of four particles: dV — d^ri2d^ri3d^ri4. 
An integral with all njk — is called "generating" : 

4 4 



i<fe 



i<k 



All the integrals in Eq. (10) can be obtained from the 
generating integral, Eq. (11), by differentiation: 



J(K,},{a,fc})-[n(- 



d 



j<k 



daj, 



r^^]I{{a,k}) . (12) 



The generating integral is given by the following formula: 



167r3 



EE-(7l^'V-) + E"(M^^M^')] 



j=i k=i 



J=2 



(13) 



The functions v{z) and u{z) are expressed in terms of the 
dilogarithmic function Li2(2;): 



u{z) = Li2{z) - Li2(l/2;) 



(14) 



v{z) = i Li2[(l - z)/2] - i Li2[(l + z)/2] - (15) 
^iln2[(l-z)/2] + iln2[(l + z)/2] . 



In Eq. (13) for the generating integral, 7^"'' are third 
order polynomials in a's, defined in the following way: 



(16) 



7. 



where for each j k: I ^ j^k] m j^k; I ^ m. The 
polynomials are defined as follows: 



(?) /2 2 2\ 

Mfc = aimi-ajk + aki + "fcm) . 
/i- = 2ahnakiak7n , 



(17) 



with the same restrictions on values of j, k, I, and m. 

The function ct is a square root of a sixth-order poly- 
nomial in a's: a — y/si + S2- The quantity si in this 
expression is given by 

4 

E 22/2,2 2 2 2 2 \ /"i o^ 

aija;m(aij+aim-ai/-aim-aj7-ajm) > (18) 

where for each j: I ^ 1, j; m ^ 1, j; I ^ m. The quantity 
S2 is determined as 



E2 2 2 



(19) 



where for each j: l^m^k ^ j] I ^ m] m ^ k] I ^ k. 

(i) 

Finally, is defined by the following expression: 



/3i" = (^-7i'^)/(- + 7i'^) 



(20) 



In all these formulas, indices j, k, I, m change from 1 to 
4, and it is assumed that ajk = otkj for each j ^ k. If 
some indices are not defined uniquely, the formulas are 
symmetric under their permutations. 

Eq. (13) is the main result of this method |Q. It pro- 
vides an analytic expression for the generating integral, 
Eq. (11). It was pointed out Q that there is no need 
to know an analytic dependence of the generating inte- 
gral on the parameters {otjk} to compute the family of 
integrals, Eq. (10). According to Eq. (12), all these in- 
tegrals are derivatives of the generating integral. Special 
formulas can be used Q to calculate numerical values 
of derivatives of functions / • g and ^1(5), if numerical 
values of derivatives of the functions /, and h have 
already been computed. For example, derivatives of the 
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term v{'^^^ /a) in Eq. (13) can be obtained in the follow- 
ing way. First, derivatives of cr^ and 7^"*^ with respect 
to {otjk} are calculated. Then derivatives of a function 



h{z) 



'1/2 



are computed at z = ct . After that, us- 



ing a formula for derivatives of h{g) with g cr , one 
finds derivatives of 1/cr with respect to {ajk\. Then, us- 
ing a formula for derivatives oi f ■ g with / = and 

g — 1/(7, one obtains derivatives of 7^"'^ with respect to 
{ctjk}- After that, derivatives of a function h{z) = v{z) 

at z — /u are calculated. Finally, using a formula 

for derivatives of h{g) with g = ^l^^ /cr, one can find the 

derivatives of v{jI^^ /a) with respect to {ajk}- Within 
this approach, all the integrals of Eq. (10) can be evalu- 
ated by means of an efficient recursive procedure, working 
with numbers only. 

At this point, we can appreciate importance of 
Eqs. (4)-(9). In order to obtain the matrix elements of 
the Hamiltonian, Eq. (1), we have to compute the mixed 
derivatives given by Eq. (12) up to the sixth order only, 
i.e. for Ujk — 0,1, where j,k — 1...4, and j < k. This 
means that, at each step of the recursive procedure, we 
calculate 2^ = 64 derivatives. If we tried to evaluate 
the integrals of Eq. (3) directly, it would be necessary 
to compute the mixed derivatives in Eq. (12) up to 18th 
order, i.e. for Ujk ~ 0...3. The number of derivatives, 
calculated at each step, would increase quadratically. An 
amount of time, required to carry out the entire recursive 
procedure, would be enormous. Therefore, the original 
method by Fromm and Hill Q], used by itself, does not 
make high precision calculations of four-particle systems 
possible. Only in conjunction with the method |12 for 
reducing the number of integrals can it produce valuable 
results. 



C. The multiple valued functions 



In the immediate vicinity of the unit circle, where con- 
vergence of the series in Eq. (22) is slow, the following 
relations can be used to shift the argument of Li2(2;): 



Li2(z) = — -ln(z)ln(l-z)-Li2(l-z) . (24) 
6 



Li2(z) = -Li2(z2)-Li2(-z) 



(25) 



Presence of the logarithm in Eqs. (23) and (24) clearly 
indicates that the function Li2(2:) is, in general, multiple 
valued. In order to specify its principal branch we need 
to fix the principal branch of the logarithm. The complex 
logarithm has branch points at and 00. We choose its 
branch cut to run along the negative real axis and define 
the principal branch as follows: 



ln(z) = In |z| + i argz , 



TT < arg z < TT 



(26) 



This choice determines branch cuts and fixes the prin- 
cipal branch for the dilogarithm, and the functions u{z) 
and v(z). 

The function Li2(z) has branch points at 1 and 00; its 
branch cut runs from 1 to cx) along the positive real axis. 
The function u{z) has branch points at 0, 1, and 00; its 
branch cut goes from to 00 along the positive real axis. 
The function v{z) has branch points at 1, -1, and 00; its 
branch cuts run from 00 to -1 along the negative real 
axis, and from 1 to 00 along the positive real axis. 

Fig. 1 exhibits the branch points and cuts for these 
multiple valued functions. 

It is important to note that the function a, which is 
present in Eq. (13) and defined using Eqs. (18) and (19), 
is also a multiple valued function. The complex square 
root has branch points at and 00. We choose its branch 
cut to run along the positive real axis and define the 
principal branch as follows: 



The main difficulty in using Eq. (13) for the generating 
integral is the fact that the functions in this formula are 
multiple valued. Indeed, the functions u{z) and v{z), 
given by Eqs. (14) and (15), are expressed in terms of the 
dilogarithmic function Li2(z). The dilogarithm is defined 
as follows [pi: 



Li2(^) 



ln(l - C) 
C 



dC 



(21) 



This function is analytic inside the unit circle in the com- 
plex plane: 



Li2(^) = E 



\z\<l 



(22) 



71=1 



Its values outside the unit circle can be determined using 
a relation [^: 

Li2(^) = -y-^ln^(-^)-Li2(l/z) . (23) 



G ^ V|z|exp(-argz) 



< argz < 27r 



(27) 



It can be seen from the definition of the generating in- 
tegral, Eq. (11), that it is a continuous function of param- 
eters {ctjk} for all values of these parameters satisfying 
the following conditions: 

ai2 + ai3 + ai4 > , Q!i2 + 0123 + "24 > , (28) 

ai3 + Oi23 + "34 > , ai4 + Q!24 + "34 > . 



"12 + "13 + "24 + "34 > 
"12 + "14 + "23 + "34 > 
"13 + "14 + "23 + "24 > 



(29) 



These conditions mean, physically, that the wavefunc- 
tion of a system of four particles decreases exponen- 
tially when any of the interparticle separations become 
infinitely large. If the parameters {ctjk} are complex, the 
above inequalities must be satisfied by their real parts. 
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a) 


b) 
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c) 






d) 
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-1 


1 



FIG. 1: Branch cuts in the complex z plane, necessary to 
define principal branches of the multiple valued functions: a) 
ln(z), branch points at and oo; b) lA^i^z), branch points at 1 
and oo; c) u(2:), branch points at 0, 1, and oo; d) f (z), branch 
points at 1, -1, and oo. 

The continuity of the generating integral, Eq. (11), im- 
plies that the right hand side of Eq. (13) is also a con- 
tinuous function of {oLjk\- This fact has two important 
consequences. 

First, the multiple valued functions u{z)^ ''^(^)i ^-nd 
<j{z) in Eq. (13) remain continuous while their branches 
change. As a point of interest moves in 12-dimensional 
space of six complex parameters, {ofjfe}, the arguments 
of these functions move freely in the complex plane, and 
their branches change repeatedly. However, a computer 
can evaluate only the principal branch of the logarithm, 
given by Eq. (26), and the principal branch of the square 
root, given by Eq. (27). Therefore, only the principal 
branches of the functions Li2(z), m(z), v(z), and (7(z), de- 
fined in the complex plane with the branch cuts, can be 
calculated directly. Thus, a special procedure of branch 
tracking is necessary to restore continuity of these func- 
tions every time their arguments cross the branch cuts. 

Second, all singularities, which different terms in 
Eq. (13) can have, cancel mutually. These singularities 
arise when tr = 0, and when any of the following equali- 
ties are satisfied: 

- OLjl -I- a^m + OLjn = , 

OL^l - ajm + Oijn = , (30) 

where for each j = 1...4: l,m,n ^ j; I ^ ni] m ^ n; 
I ^ n. These singularities are unphysical, and should 
have no effect on the value of the generating integral. 
As a point under consideration moves in the space of 
the parameters {aijk}^ the arguments of the functions 



u{z) and v{z) can frequently appear in the vicinity of 
the singular (branch) points. As a result, the values of 
these functions can exhibit considerable change, even if 
the parameters {(Xjk} change only slightly. Therefore, 
a special procedure for dealing with the singularities is 
needed in order to carry out explicit cancellation of all 
singular terms. 

This discussion demonstrates that the method of Q| 
is impossible to use without an effective algorithm for 
numerical branch and singularity tracking. 

D. Branch tracking in the real case 

Before discussing a general algorithm of branch track- 
ing, it is beneficial to consider a particular case, when 
all the exponential parameters, {a^fc}, are real numbers. 
Let us introduce the following parametrization: 

ajkip) = (ajk -l)p+l , < p < 1 . (31) 

As the real parameter p changes from to 1, the cor- 
responding point in 6-dimensional space moves from 
(1,1,1,1,1,1) to (ai2,ai3,ai4,a23,a24,Q;34)- If the pa- 
rameters {ctjk} satisfy the conditions of Eqs. (28) and 
(29), the parameters {ajk{p)} will satisfy these condi- 
tions for any p between and 1. Therefore, the gen- 
erating integral, given by Eq. (13), must be a con- 
tinuous function of p. It is known [Q that Eq. (13) 
with the functions u{z), v{z), and cr(z), represented by 
their principal branches, yields the correct value for the 
generating integral at the reference point (1,1,1,1,1,1). 
If this value changes continuously, as the parameter p 
goes from to 1, one can be sure that the generating 
integral will be computed correctly at the final point 
(ai2, aia, 0:14, 0:23, 024, 034). Therefore, continuity of this 
integral is a criterion of the correct branch tracking. 

Let us define a function S{p) as the sum in the square 
brackets of Eq. (13) when the parameters {ajk{p)} are 
used instead of {ctjk}- 

^(?')-EE-(7l''V-)+E«(M'^/3?'^) . (32) 

J = l k=l 3=2 

Our purpose is to ensure that this function is continuous 
along the path from p — to p = 1. 

First, we consider a case when ct^(p) > 0. The func- 
tion a{p) is real, and all the arguments of the functions 
u{z) and v{z) in Eq. (32) are real as well. It will be 
shown in Sec. ILE that only imaginary parts of these 
functions exhibit discontinuities, when their arguments 
cross the branch cuts. Because the generating integral is 
real, the imaginary parts of the functions u{z) and v{z) in 
Eq. (32) must cancel anyway. Therefore, discontinuities 
in the real part of S{j)) may appear near the singular 
points of the functions u{z) and v{z) only. The singu- 
larities of different terms in Eq. (32) should cancel one 
another. However, because of possible branch changes. 
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complete cancellation may not happen. The formulas of 
Sec. II. E suggest that, near the singular points of u{z) 
and v{z), the real part of the function S{p) can undergo 
changes by m7r^, where m is some integer. Thus, the 
function S{p) can have finite discontinuities, which are 
integer multiples of tt^. 

From now on, the branch tracking is only a techni- 
cal problem. To solve it, it is necessary to find all val- 
ues of the parameter p between and 1, which cor- 
respond to singular points. They include zeros of the 
sixth-order polynomial (p) , and values of p, at which 
the parameters {ajk{p)} satisfy any of the conditions of 
Eq. (30). Let us denote the resulting set of numbers as 
{Pjli j — I------ The correction function, needed to re- 
move discontinuities of the function S{p), is given by the 
following expression: 

C{p) - -TT^ Nint[(5(p, +e) - S{p, - e))/^^] _ (33) 
Pj<P 

Here, the function Nint[x] returns an integer number, 
nearest to the real number x. The value of e in actual 
calculations was set to 10~^. The correct value of the 
generating integral can now be determined from the for- 
mula: 

^ j=l k=l j=2 

(34) 

Therefore, in the case of the real parameters {ctjk}, 
the procedure of branch tracking can be implemented 
without a detailed numerical analysis of behavior of the 
multiple valued functions. All we need to do is to calcu- 
late the function S{p) twice for each singular point pj, 
encountered along the path from p — to p = 1, and 
subtract discontinuities, proportional to tt^. The time, 
needed to determine the correction C(l) in Eq. (34), is 
shorter than the time, required to carry out the recursive 
procedure for the family of integrals. It does not increase 
the overall computation time significantly. 

The case of cr^(p) < is also straightforward. The 
quantity a is now imaginary. The function S{p) is imag- 
inary as well, thus giving a real value of the generating 
integral. Im(5(p)) can be expressed in terms of Clausen's 
function Cl2(0), which is a real function of a real argu- 
ment Q 1^. Eqs. (33) and (34) are valid also in this case, 
if S{p) is replaced by lm{S{p)), and a{p) is replaced by 
lni{a{p)). Therefore, in both cases (cr^ > and < 0) 
the entire algorithm for analytic evaluation of the four- 
particle integrals can be presented in the real form with- 
out any use of complex numbers. 

This simplified method of branch tracking has been 
successfully employed in variational calculations of four- 
particle systems |^. Therefore, the described method 
of branch tracking in the case of real exponential param- 
eters is both theoretically correct and practically reliable. 



E. Branch tracking in the complex case 

Let us now describe a method of branch tracking in 
a general case, when the exponential parameters, {a^fe}, 
are complex numbers. It is assumed that their real parts 
satisfy Eqs. (28) and (29). We use the same parametriza- 
tion as before, but with a complex parameter p : 

«jfe(p) = (a^-fe - 1)P + 1 , < Rc{p) < 1 . (35) 

As p moves in the complex plane from to 1, 
the corresponding point in 12-dimensional space of 
six complex parameters moves from (1,1,1,1,1,1) to 
(ai2, ai3, ai4, a23, ^24, 034)- The generating integral, 
Eq. (13), must be a continuous function oi p. Moreover, 
its value, computed at the final point, {ajk}, should not 
depend on a choice of the path from p = to p = I. How- 
ever, an optimal choice of this path can facilitate branch 
tracking considerably. 

Fig. 2 exhibits three examples of paths in the complex 
p plane. In case a), there are no singular points on or 
near the real axis between and 1. The path is simply 
a straight lime segment between these points. In case 
b), there is one point, pi, at which different terms in 
Eq. (13) exhibit singular behavior. The path is the same 
as before, except for a small semicircle in the vicinity of 
this point. In case c), there are two singular points, pi 
and p2, near the real axis between and 1. The path 
is more complicated, as shown in the figure. In general, 
only those singular points in the p plane, which are close 
to the real interval between and 1, are of interest. The 
path should be carefully defined in the vicinity of every 
such point to allow precise analysis of behavior of all 
arguments of the multiple valued functions. The values of 
p, at which singularities may arise, can be found from the 
polynomial equation <J^{p) = 0, and from twelve linear 
equations, contained in Eq. (30). 

In order to obtain correction functions for the function 
u{z), defined by Eq. (14), we have to consider behavior 
of this function near its branch points 0, 1, and 00: 

1 2 

u{z^O)= -In (-z) + U(o)(2;) ; 

m(z ^ 1) = -21n(z)ln(l - z) + U(i)(z) ; (36) 
1 2 

U{Z ^ 00) = -- In (-Z) + M(oc.)(z) . 

In these formulas, the functions with subscripts (0), (1), 
and (00) are functions, analytic in the vicinities of 0, 1, 
and 00, respectively. 

Let us introduce the following notations. A com- 
plex function z{p) will represent any of the arguments, 

of the fimction u(z) in Eq. (13). It depends 
on p through the parameters ajk(jp), given by Eq. (35). 
Let {pj}, j — 1...N, denote values of the parame- 
ter p, for which z{pj) are singular points 0, 1, 00, or 
any points, where z{p) crosses the real axis. It is as- 
sumed that < Ke{pj) < 1 for each j = 1...N, and 
Re(pj) < Re(pj+i). Each point pj will be characterized 
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If z{p) moves near the singular point z{pj) =00, 
let rij =5, and 

U5{z,j) = +5|/2 + i Sj[\n{-z) + Uj] . 



(41) 



If z{p) crosses the branch cut ] — oo,0 [ of the 
function ln(z) at z{pj), let rij = 6. 

The logarithms in these formulas are multiple valued 
functions themselves. Their branches can also change, 
and they can exhibit singular behavior, while an argu- 
ment z{p) moves further in the complex plane. Because 
only the principle branch of the logarithm is calculated 
by a computer, the additional terms, Uj and Uj are in- 
cluded to correct values of these functions. These terms 
are given by the following formulas: 



JV 



N 



Uj = + ^ 2mfe7r i - ^ iSk 



(42) 



nfc = l,2 



k>j 
rafc=4,5 



FIG. 2: Examples of paths in the complex p plane: a) no 
singular points on or near the real axis between and 1; b) 
one singular point, p 1, on the real axis; c) two singular points, 
pi and p2, near the real axis. The plots are not to scale. 



by an index nj, and either integer nij, or real Sj. The 
index nj ~ 1...6 specifies a type of singular behavior, 
as explained below. The number rrij provides informa- 
tion about direction, in which the real axis is crossed by 

z(p). We set rrij = -1-1, if the axis is crossed from be- 
low (i.e. t)j Emd rrij ~ —1, if it is crossed from above 
(i.e. I). The real quantity Sj is equal to a change in 
arg(z(p) — z{pj)), when z{p) moves in the vicinity of a 



singular point z{pj). If z(j)j) 



the quantity 6j de- 



notes a change in a,rg{z{p)). These notations will allow 
us to present the algorithm of branch tracking as a series 
of formulas. 

Five correction fimctions, Un-{z,i), are needed to re- 
store continuity of the computed function u{z). 

If z{p) crosses the branch cut ] 1, -|-oo [ at z{pj), 
let rij = 1, and 

ui{z,j) = +27r^ - 2mj7r i [ \n{-z) + Uj] . (37) 

If z{p) crosses the branch cut ] 0, 1 [ at z(pj), 
let rij =2, and 

U2{z,j) = -27^"^ + 2mjTTi[hi{-z) + Uj] . (38) 

If z{p) moves near the singular point z{pj) = 1, 
let rij =3, and 

u^{z,j) = 2i5j[\n{z) + JJj] . (39) 

If z{p) moves near the singular point z{pj) = 0, 
let rij =4, and 



JV 



N 



(43) 



k>j 
nk=6 



k>j 
nji=4,5 



The condition rifc = 1 , 2 in these formulas means that we 
have to sum up only those indices ruk, that correspond 
to situations, when z{p) crosses the branch cuts ] 1, -l-oo [ 
and ] 0, 1 [. The condition n/j = 4, 5 limits the summation 
of 5k to those cases, when z{p) moves near the singular 
points and 00. If rik = 6, we consider only situations, 
when z{p) crosses the real axes in the interval ] — 00, [. 

Thus, each singular or crossing point z{pj), encoun- 
tered by the argument z{p) of the function u{z), gives 
rise to a correction fimction, Un {z, j), required to make 
u{z) continuous. However, the structure of this correc- 
tion function at the end of the path, p = 1, will depend 
on behavior of z{p) near all the following singular and 
crossing points, z{pk), j < k < N . The resulting correc- 
tion function, Uc{z), obtained after passing all the points 
z{Pj), j = 1---N, is given by the following expression: 



N 



(44) 



In order to see, how these correction functions operate, 
consider values of the principal branch of u{z) at the 
edges of the branch cut ] 1, -|-oo [: 

u{x±ie) = —-2U2{l/x)--ln^{x)±iTrln{x) . (45) 

o ^ 

In this formula, a: > 1 is real, and e +0. Imagine that 
the branch cut is crossed from below (t)- Then rrij = -1-1, 
and the value of the correction function ui{z,j), defined 
by Eq. (37) , at the point z = x + ie is equal to 



U4{z,j) = -6^/2 - i6j[ln{-z) + Uj] 



(40) 



ui{x + ie, j) = —2Tr i ln{x) . 
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If the branch cut is crossed from above (J,), rrij = —1. 
and the value of this correction function at the point 
z = X — ie\s equal to 

u\{x — i€,j) = +2-Ki ln(x) . 

Thus, the correction function ui{z, j), added after the 
branch cut is crossed, eliminates the finite discontinuity 
of the principal branch, Eq. (45), of the function u{z) 
along ] 1, +00 [. The correction function U2{z,j), defined 
by Eq. (38), acts in a similar way at ] 0, 1 [. 

Imagine now that the argument of u{z) goes around 
the singular point at oo, starting from z = x + ie, and 
coming back to z = x — ie, without crossing the branch 
cut along the positive real axis. The value of u{z) exhibits 
a singular change (from Eq. (45)) by 

Au = — 27ri ln(a;) . 

In this case, 6j = 27r, and the value of the correction 
function 1*5(2, j), defined by Eq. (41), at the point z = 
X — ie is equal to 

U5{x — ie,j) = +2'7ri ln{x) . 

Thus, the correction hmction U5(z,j), added after z has 
moved near the singular point at 00, eliminates the sin- 
gular contribution to the value of the function u{z). 
The correction functions u:i{z, j) and U4(z, j), given by 
Eqs. (39) and (40), produce similar results for the other 
singular points. 

If, in the above examples, the argument of u{z) first 
crosses the branch cut, and then moves around the 
singular point, the correction function ui{z,j) has to 
be modified by adding nonzero Uj to the logarithm 
according to Eq. (37). 

The same principles of branch tracking apply to the 
function ^^(2;), defined by Eq. (15). First, we consider 
behavior of this function near its branch points 1, -1, 
and 00: 

viz^ ^) = -\^Ai^-z)/{l + z)]+v^,^{z) ; 

v{z^-l)= iln2[(l + z)/(l-z)]+t;(_i)(z) ; (46) 

v{z ^ 00) = i ln(-zV4) ln[{z + l)/{z - 1)] + v^^^{z) . 

In these formulas, the functions with subscripts (1), (-1), 
and (00) are functions, analytic in the vicinities of 1, -1, 
and 00, respectively. 

Let us again consider a complex function z{p), which 

can represent each of the arguments j^^^ /a of the func- 
tion v{z) in Eq. (13). Let {pj}, j = 1...N, denote values 
of the parameter p, such that z{pj) are singular points 1, 
-1, 00, or z{p) crosses the real axis. It is assumed that 
their real parts form an increasing set of numbers be- 
tween and 1. As before, each point pj is characterized 
by an index rij , and either mj or 6j . Values of rij will be 



assigned below, and meanings of mj and Sj remain the 
same. 

Five correction functions, Vnj{z,j), are used to make 
the computed function v{z) continuous. 

If z{p) crosses the branch cut ] 1, +00 [ at z{pj), 
let rij = 1, and 

v^{z,j)=W+mjni{\u[{l^z)/{l-z)] + Vj} . (47) 

If z{p) crosses the branch cut ] — 00, — 1 [ at 2;(pj), 
let rij =2, and 

V2{z,j) = -TT^-mjTTi{\n[{l + z)/{l-z)] + Vj} . (48) 

If z{p) moves near the singular point z{pj) = 1, 
let Uj =3, and 

vs{z,j) = +1 - i |{ ln[(l + z)/{l - z)] + Vj}. (49) 

If z{p) moves near the singular point z{pj) = —1, 
let rij =4, and 

v,{z,j) = -| - i |{ ln[(l + z)/{l - z)] +Vj}. (50) 

If z{p) moves near the singular point z{pj) =00, 
let rij =5, and 

vr,{z,j) = -i5j{\n[{z + l)/{z-l)\ + Vj} . (51) 

If z{p) crosses the branch cut ] — 1, 1 [ of the 
function \d\{z + l)/{z — 1)] at z{pj), let rij = 6. 

The additional terms, Vj and Vj, necessary to correct 
behavior of the logarithms, are given by the following 
formulas: 



N 


N 


N 


2mkTTi ^ 




E 


k>j 


k>j 


k>j 




nk=3 


nk=i 


N 


N 


N 



t^- = + E 2mfe7ri + ^ i4- E ^^'^ ■ (^"^) 

k>j k>j k>j 

nfc=6 nk=3 nfc=4 

As in the previous case, a correction function Vnj{z,j) 
has to be added to the function v{z) every time its argu- 
ment z{p) passes a singular or crossing point z{pj). This 
way, the calculated function v{z) can be made continu- 
ous. However, the form of this correction function at the 
end of the path depends on behavior of z{p) near all the 
points z{pk), following z{pj). The resulting correction 
function, Vc{z), is the following: 

JV 

Mz) = ^Vnj{z,j) . (54) 



9 



Let us now briefly discuss the effect of using these cor- 
rection functions. Consider values of v{z) at the edges of 
the branch cut ] 1, +00 [: 

v{x ± ie) = iLi2[2/(l + x)] - iLi2[2/(l - x)] + 

+ iln^[2/(l + a;)]-iln2[2/(a;-l)]± (55) 

±i^H{x-l)/{x + l)] . 

Here, a; > 1 is real, and e +0. Imagine that the branch 
cut is crossed from below ("f ). Then mj = +1, and a value 
of the correction function, Vi{z,j), defined by Eq. (47), 
at a point z = x + ie is equal to 

Vi{x + ie,j) = — « 7rln[(.T — 1) / {x + 1)] . 

If the branch cut is crossed from above (|), then mj = 
— 1, and a value of this correction function at z = x — ie 
is equal to 

v-i{x — ie,j) = +in\n[{x — l)/{x + 1)] . 

Therefore, the function vi{z,j), added to the function 
v{z) after the branch cut is crossed, removes the discon- 
tinuity of the principal branch along ] 1, -|-oo [. The cor- 
rection function V2{z,j), given by Eq. (48), makes v{z) 
continuous at ] — 00, — 1 [. 

Imagine now that the argument of v{z) moves around 
the singular point +1, starting from z = x + ie and re- 
turning to z = X — ie, without crossing the branch cut. 
The value of v{z) undergoes a change (from Eq. (55)) by 

Av = - i7rln[(a; - l)/(x -(- 1)] . 

Because 6j = 2tt, a. value of the correction function 
vz{z,j), defined by Eq. (49), &t z = x — ieis equal to 

V'ilx — ie,j) — +i'Kh\[{x — l)/{x + 1)] . 

Thus, by adding the correction function V3{z,j), it is pos- 
sible to eliminate the singular contribution to the value of 
v{z), when z goes around the singular point +1. The cor- 
rection functions VA{z,j) and v^{z,j), given by Eqs. (50) 
and (51), produce the same results for the other two sin- 
gular points. 

If, in the above examples, the argument of v{z) first 
crosses the branch cut, and then moves around the 
singular point, the correction function i>i(z, j) should be 
modified by adding nonzero Vj according to Eq. (47). 

It is important to note that the function a{z) is also a 
multiple valued function. Its principal branch, defined by 

Eq. (27), changes sign each time the argument z crosses 
the branch cut along the positive real axis. If this hap- 
pens N times while the parameter p changes from to 
1, the corrected value (Jc{z) of this function at p = 1 is 
equal to 

a,{z) = {-lfa{z) , (56) 



where a{z) is the value of the principal branch of the 
complex square root. 

We are now in a position to write a corrected expres- 
sion for the generating integral: 

I = ^[EE(-+-c)(7i^'V-c)+i:(-+-c)(/3r^/3F^)] 

j=l k=l 3=2 

(57) 

In this formula, each of 16 terms with the function v[z) 
contains its own correction function fc(z), which is char- 
acterized by its own and sets of numbers {pi}, {ni}, 
{mi}, and {5i}. The same is true for each of the three 
terms with the function u{z). 

Eq. (57) is profoundly different from the original for- 
mula, Eq. (13), for the generating integral. In Eq. (13), 
the functions u{z), v{z), and cf{z) are expressed in terms 
of the multiple valued logarithm and square root. When 
Eq. (57) is used, it is assumed, on the contrary, that 
all the logarithms and square roots are represented by 
their principal branches, and, therefore, can be readily 
evaluated by a computer. The multiple valued nature 
of the functions u{z), v{z), and a{z) is taken into ac- 
count explicitly by means of the additional correction 
terms and factors. Also, singular contributions from 
different terms in Eq. (13) are expected to cancel each 
other to yield a correct value for the generating integral, 
Eq. (11). When we use Eq. (57), all the singular con- 
tributions from different terms are cancelled explicitly 
and separately, so that each hmction [v + Vc){'^\^^ /(^c) 
or [u + Uc)(/3i^^/3i''^) is continuous along the path from 
(1,1,1,1,1,1) to (a;i2,ai3,ai4,a23,Q;24,Q!34). As a result, 
the generating integral, obtained from Eq. (57), is a con- 
tinuous function of the complex parameters {cc^fe}. Thus, 
the problem of branch tracking is successfully solved in 
the most general case. 

Let us consider an example of branch tracking for the 
function v{z) with the argument z{p). A sample path 

in the complex p plane is displayed in Fig. 3a, and the 
corresponding path from z(0) = A to z(l) = B in the 
complex z plane is shown in Fig. 3b. There are five points 

of interest: z{pi) = 1, ^(p;?) = 00, and z{py,) = —1 are 
singular points of the function v{z): z(p2) and z{p4) are 
points, where the argument z{p) crosses the branch cuts. 
According to the chosen classification: ni =3, i5i = tt; 
n2 = 1, m2 = +1; = 5, ^3 = tt; 714 = 2, 1714 = — 1; 
ns = 4, 65 = n. Using Eqs. (47)-(53), one can easily 
obtain expressions for the correction functions: 

2 

vsiz,!)^ z|{ln[(l + z)/(l-z)]- 7 7r} ; 

vi{z,2)= + iTr{\n[{l + z)/{l- z)]+ in} ; 
V5{z,3)= - in{\n[{z + l)/{z-l)]- in} ; 
V2{z, 4) = -n'^+ i n{ ln[(l + z)/{l -z)]- in} ; 

V4{z,5) = -^-i^{H{l + z)/{l-z)]+ 0} . 




FIG. 3: Illustration of branch tracking in the complex case: 
a) the path in the complex p plane, chosen to avoid singu- 
lar points, pi, p3, and ps; b) the corresponding path in the 
complex z plane from z(Q) =Ato z{l) = B, where z{p) is an 
argument of the multiple valued function v{z). 

The resulting correction function, Vc{z), is 

Vc{z) = ^+ i7rln[(l + z)/(l-2:)]- 

-i7rln[(z + l)/(z- 1)] . 

One can see that this function is different from zero, even 
if yl = _B, i.e. the contour is closed. This is not surpris- 
ing. Even though v{z) is represented by its principal 
branch, the sum v{z) -\- Vc{z) is still a multiple valued 
function. Its value generally undergoes a finite change, 
if its argument traverses a closed loop, encircling branch 
points. Consider a value of Vc{z) a,t z — x — ie, where 
— 1 < a; < 1, and e -1-0. In this case, the logarithms 
cancel, and 

Vc{x — ie) = — 7r^/2 . 

Contributions of this type from different terms in Eq. (57) 
produce an additional constant ttitt^, needed to correct 
a value of the generating integral in the case of real pa- 
rameters {otjk\- Thus, Eq. (34) is a particular case of 
Eq. (57). 

This example demonstrates that the branch tracking 
in the general case requires a comprehensive numerical 
analysis of behavior of all the arguments in Eq. (13). 

F. Numerical procedure and results 

Practical implementation of the method, described in 
the previous sections, is inevitably a very complicated 




FIG. 4: Segmentation of the path in the complex p plane, 
needed to analyze behavior of arguments of the multiple val- 
ued functions numerically. 



task. Detailed information about the recursive proce- 
dure, needed to compute the family of integrals, Eq. (10), 
can be found in Q . Here we describe only the procedure 
for numerical branch tracking. 

First, the set of points, {pfe}, at which different terms 
in Eq. (13) can exhibit singular behavior, is determined. 
This is done by solving the sixth-order equation ct^ = 0, 
and linear equations of Eq. (30), with the parametriza- 
tion according to Eq. (35). Only those values of p, that 
lie in or near the real interval ] 0, 1 [, are included in the 
set {pk\- Then, a path from to 1 in the complex p 
plane is chosen. Fig. 2 gives an idea of this. The whole 
path is shifted downward by a small imaginary quan- 
tity ie to avoid possible ambiguities, when arguments of 
the functions u{z) and v{z) are real. All the arguments 
in Eq. (57) are computed at the final point of this path, 
p = 1 — ie. In actual calculations, e was set to 10~^®. This 
did not affect values of the integrals, but was enough to 
shift the arguments from the real axis. 

The path in the complex p plane is divided into small 
intervals, as shown in Fig. 4. The intervals ]P/,P;+i [, 
into which the linear segments between the singular 
points, {pfe}, are divided, have a typical length of 10~^. 
Each small semicircle beneath a singular point, pk, has 
a radius r = 10^^, and divided into six parts. The cor- 
responding boundary points are 

-P; =Pfc +7-exp[i (7rZ/6)] - ie , Z = 0...6. (58) 

In order to obtain full information about behavior 
of different arguments in Eq. (13), all these arguments 
should be computed at all the points Pi along the path. 
The quantities 7^-'^ (3\^\ and ct^, given by Eqs. (16)-(20), 
are simple functions of {ajk{p)}^ so this calculation can 
be performed almost immediately. The values of each 
argument are analyzed, and the numbers N , nj, rrij, and 
5j, j = l...iV, needed to apply the formulas of Sec. II. E, 
are determined. This procedure works as follows. To find 
out, if an argument crosses the real axis, the imaginary 
parts of its values, computed at points Pi and are 
compared. If they have opposite signs, dichotomy is used 
to reduce the interval, and determine, where the real axis 
is crossed, and in which direction. This is also done for 
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the intervals along each small semicircle, but without the 
dichotomy. In this manner, all the crossing points can be 
found and analyzed. At each interval, ]Pi,Pi+i [, this 
analysis is carried out for cr^ first. If cr^ crosses the posi- 
tive real axis, when p is between Pi and Pi+i, the sign of 
a should be changed. Then, all computed at Pi+i, 

should be multiplied by (-1), and all - inverted. The 
quantity Sj is determined as follows: 

5 

<5,=^[arg(2(P,+i))-arg(z(P,))] • (59) 
/=o 

Here, the points Pi are specified by Eq. (58), and z can 
stand for any of the arguments "fi^^ /cr and as 
functions oi {ajk{p)}- The right hand side of this formula 
is presented as a sum, because as p changes from Pq to 
Pe, the argument z{p) may go around the singular point, 
z{pk) several times. Thus, behavior of each argument 
in the vicinity of each singular point can be analyzed. 
The described procedure of numerical branch tracking 
provides all the information, necessary for successful use 
of Eq. (57). 

The entire algorithm for analytic evaluation of the 
four-particle integrals with complex parameters was 
tested in four different ways. 

First, real parameters {ajk} were used, and the re- 
sults, obtained using the complex algorithm of Sec. II. E, 
were compared with results, provided by the method of 
Sec. II. D for the real case. The real parts of the computed 
integrals, Eq. (10), were invariably in excellent agree- 
ment. The imaginary parts, given by the new method, 
were at least 20 orders of magnitude smaller than the 
real parts, and could be considered negligible. Therefore, 
the complex algorithm works correctly for any acceptable 
real parameters. 

Second, the parameters {ctjk} were multiplied by an 
arbitrarily chosen complex number A. A resulting inte- 
gral, Eq. (10), with a particular set {rijk} must be equal 
to its original value, multiplied by / = — 1/(— A)^, where 
K = ni2+ni3 + ni4+n23 + n24:+n34 + 3. Various A's were 
used, and values of the integrals, calculated directly, were 
compared with the rescaled original values. Remarkable 
agreement was observed in all these cases. Note that dif- 
ferent values of A correspond to different paths in the 
space of parameters according to Eq. (35). 

Third, if two exponential parameters, ais and 024, are 
equal to zero, the six Coulomb integrals and one over- 
lap integral, needed to determine matrix elements of the 
Hamiltonian according to Eqs. (4)- (9), can be obtained 
analytically in terms of rational functions and logarithms. 
Values of these integrals, calculated with various sets of 
complex parameters, ai2, ai4, 023, 0:34, were compared 
with the same integrals, computed using the new method. 
They were always in complete agreement. 

Fourth, different paths in the complex p plane were 
chosen. They included singularities, located not only 
near the real axis, but also further away. The results did 



not depend on the choice of the path. This fact suggests 
that the described method of numerical branch tracking 
is stable and reliable. Of course, the path in actual cal- 
culations should be as simple as possible, provided that 
all nearby singularities are carefully taken into account. 

Table I displays values of the integrals for three differ- 
ent sets of parameters {ajk}, used to test the computer 
program. Many other sets of parameters were also con- 
sidered. All the integrals were calculated using the gen- 
eral algorithm for numerical branch tracking, described 
in Sec. II. E. Implementation of this algorithm requires 
quadruple precision. The program computes a family of 
64 integrals, Eq. (10), with two possible values for ev- 
ery index: rijk = 0, 1. Only seven integrals, necessary 
to obtain matrix elements of the Hamiltonian according 
to Eqs. (4)-(9), are presented in Table I for each set of 
parameters. 

Our results demonstrate that the developed algorithm 
allows precise evaluation of the four-particle integrals 
with arbitrary complex parameters, provided that the 
integrals themselves converge. 

The described method makes it possible to use the 
highly versatile exponential-trigonometric basis functions 
in variational calculations of four-particle Coulomb sys- 
tems. In order to illustrate efficiency of the new basis, we 
would like to mention some results, obtained previously 
|ll| for the following systems: e'^e~e'^e~, p'^ fj,~ p^ , 
/t^ e^/i+e^, and p+e~p"*"e~. The calculations were per- 
formed using one exponential-trigonometric basis func- 
tion: 

4 4 
^ Sexp{~^Ajkrjk)siii{J2Bjkrjk + C) . (60) 
j<k ]<k 

This function includes 12 real nonlinear parameters, 
{Aj-fe} and {i?j7c}, and one linear parameter, tan(C). 
It can be considered a linear combination of two expo- 
nential functions, Eq. (2), with the complex parameters 
Ajk ± i Bjk ■ The operator S ensures that this function 
has correct symmetry with respect to permutations of 
particles. 

All integrals, necessary to determine matrix elements 
of the Hamiltonian, Eq. (1), with the function ^f, were 
computed according to the method, described in this pa- 
per. The nonlinear parameters were subjected to careful 
gradient optimization. For more details about this cal- 
culation, see [pl| . 

Table II exhibits values of the ground-state energy, E, 
for e~^e~e~^e~ , p^fi~p'^^~ , /i+e"^+e~, and p+e~p"'"e~, 
determined using the variational method with the trial 
function The table also displays the most accurate en- 
ergy values, Eq, available for these systems |jl|, ^, p^ . 
One can see that the relative errors are 0.2%, 0.7%, 
2.4%, and 3.6%, respectively. The results for two adi- 
abatic systems, ^"'"e~/i"'"e~ and p~^e~p~^e~, with very 
low mass ratios, m/M, are very impressive. Neither 
Gaussian, nor exponential functions are even nearly as 
efficient [O. Thus, a single symmetrized exponential- 



12 



TABLE I: Examples of four-particle integrals evaluated using the algorithm for numerical branch tracking in the complex case. 











{n^k} 


Re (J) 


Im(J) 


ai2 




1.56 




011111 


0.20550889174003868108D+03 


-0.52323042463487687803D - 21 


Ql3 




-0.69 




101111 


0.49701602825033406834D+02 


0.37462639365280834442D - 21 


Ql4 




2.71 




110111 


0.33278606420131925558D+03 


0.29262335382214159490D - 21 


a23 




1.75 




111011 


0.69156207020089168507D+02 


0.39044885805975934492D - 21 


Q24 




1.42 




111101 


0.20401147280211976836D+03 


0.44225408142944401598D - 21 


034 




-0.50 




111110 


0.50046583501959809463D+02 


0.75995078405826138804D - 21 










mill 


0.21644781505854395857D+03 


0.32039794574654543500D - 20 


Q12 




1.56 


* (1 + Q.5i) 


011111 


-0.70977575942226172269D+02 


0.45253255249692588012D+02 




= 


-0.69 


* (1 + 0.5i) 


101111 


-0.17165677159247121876D+02 


0.10944340655611068217D+02 


cei2 


— 


2.71 


* (1 + 0.5 i) 


110111 


-0.11493589374343266153D+03 


0.73279810811752133344D+02 


ai2 




1.75 


* (1 + 0.5i) 


111011 


-0. 23884805635825330948D+02 


0.15228263175782374191D+02 


ai2 




1.42 


* (1 + 0.5 i) 


111101 


-0.70460405295819730405D+02 


0.44923522162040663029D+02 


ai2 




-0.50 


* (1 + 0.5i) 


111110 


-0. 17284824763945988644D+02 


0.11020305731851711925D+02 










mill 


-0.40739676150047588287D+02 


0.68031854740817630022D+02 


cei2 




1.29 + 


1.19i 


011111 


0.41299141847575234393D+01 


-0.13354699318829522025D+01 











101111 


0.25568585205838373519D+01 


-0. 14420903111 112389762D+01 


ai4 




2.53 - 


1.32i 


110111 


0.39327459787814931363D+01 


-0.64522069912295906967D+01 


Q23 




1.86 + 1.44i 


111011 


0.56761553854820761165D+01 


-0.12198339708832004631D+01 


Q24 









111101 


0.38294186820466745046D+01 


-0.25317810705579310712D+01 


a34 




0.65 - 


0.93 i 


111110 


0.26044120278339787630D+01 


-0.23151954768733032766D+01 










mm 


0.37849264531713841033D+01 


-0.28372153117607227596D+01 



TABLE IL Ground-state energy, E, of four molecules, com- 
puted with a single exponential-trigonometric basis function. 
The most accurate values, Eo, of this energy are taken from 
Refs. 1, 2, 14 and 15, respectively. Atomic units. 



System 


m/M 


E 


Eo 


Error 




1 


-0.514956 


-0.516003 


0.2% 


p+fi-p+fj.- 


0.1126095 


-198.2056 


-199.6294 


0.7% 




0.0048363 


-1.113198 


-1.141000 


2.4% 


p'^e~p'^e~ 


0.0005446 


-1.122378 


-1.164025 


3.6% 



trigonometric basis function, Eq. (60), provides a remark- 
able accuracy in variational calculations of various four- 
particle systems. 



III. CONCLUSION 

The method for analytic evaluation of four-particle in- 
tegrals with complex parameters, described in this paper, 
can be regarded as both further theoretical development 
and practical implementation of the original method by 
Fromm and Hill Q. Validity of this method is not lim- 
ited to the case of real parameters. Moreover, because 
the integrals are expressed in terms of multiple valued 
complex functions, it is more natural to consider a gen- 
eral case when all the parameters are complex. The orig- 



inal formula, Eq. (13), for the generating integral can be 
used only in the immediate vicinity of the standard refer- 
ence point where all the parameters are equal to 1. The 
procedure of numerical branch tracking, proposed in this 
paper, allows computation of the integrals at any other 
point in the space of six complex parameters, by taking 
into account all branch changes along the path. The sim- 
plified method of branch tracking for real parameters is 
also discussed. 

The new method makes possible high precision varia- 
tional solution of the Coulomb four-body problem in the 
basis of exponential-trigonometric functions. The first 
calculations have shown high efficiency of this basis . 
They have also demonstrated correctness of the branch 
tracking algorithm, described in this paper. 

However, if the full potential of the exponential- 
trigonometric basis is to be revealed, an efficient 
procedure for selecting optimal values of the nonlinear 
parameters is necessary. Ideally, all the parameters 
should be chosen a priori, and all matrix elements 
- computed only once. Such a procedure has been 
developed by the authors for the case of adiabatic 
three-particle systems [ p^ . All nonlinear parameters 
of the exponential-trigonometric functions had been 
chosen before the computation, which yielded 10 correct 
significant figures for the ground state energy of H2 [0 . 
We believe that the exponential-trigonometric basis can 
provide similar precision in calculations of four-particle 
systems. 
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